#include <kuai/misc/error.hpp>
#include <kuai/mol/mol.hpp>
#include <kuai/mol/molid.hpp>
#include <kuai/mol/molio.hpp>


namespace kuai {

	bool MDLMolReader::update(std::istream& stream, PhaseSpace& result) const {
		String line;
		getline(stream, line);
		getline(stream, line);
		getline(stream, line);
		if (!getline(stream, line) || line.size() < 6) {
			return false;
		}

		Index nAtoms = lexical_cast<Index>(strip(line.substr(0, 3)));
		Index nBonds = lexical_cast<Index>(strip(line.substr(3, 3)));

		std::vector<Index> ct0;
		ct0.reserve(2+nAtoms + 3*nBonds);
		ct0.push_back(nAtoms);
		ct0.push_back(nBonds);

		std::vector<XYZ>   coords(nAtoms);
		std::vector<Index> iso(nAtoms, 0);
		std::vector<Index> rad(nAtoms, 0);
		std::vector<Integer> charge(nAtoms, 0);
		std::vector<Index> stereoA(nAtoms, 0);
		std::vector<Index> implicitH(nAtoms, 0);
		std::vector<Index> valence(nAtoms, 0);
		std::vector<Index> aam(nAtoms, 0);
		std::vector<Index> inversion(nAtoms, 0);

		std::vector<Index> order(nBonds, 0);
		std::vector<Index> stereoB(nBonds, 0);
		std::vector<Index> center(nBonds, 0);

		for (size_t i = 0; i < nAtoms; ++i) {
			if (getline(stream, line)) {
				coords[i].x = lexical_cast<RealNumber>(strip(line.substr(0, 10)));
				coords[i].y = lexical_cast<RealNumber>(strip(line.substr(10, 10)));
				coords[i].z = lexical_cast<RealNumber>(strip(line.substr(20, 10)));
				String symbol = line.substr(31, 3);
				boost::algorithm::trim(symbol);
				ct0.push_back(Element::getElement(symbol.c_str()).number);
				if (line.size() >=35) {
					// TODO:
					iso[i] = lexical_cast<Integer>(strip(line.substr(34, 2)));
				}
				if (line.size() >=39) {
					charge[i] = lexical_cast<Integer>(strip(line.substr(36, 3)));
					if (charge[i] == 4) {
						rad[i] = 2;
					}
					else if (charge[i] != 0) {
						charge[i] = 4 - charge[i];
					}
				}
				if (line.size() >= 42) {
					stereoA[i] = lexical_cast<Index>(strip(line.substr(39, 3)));
				}
				if (line.size() >= 45) {
					implicitH[i] = lexical_cast<Index>(strip(line.substr(42, 3)));
				}
				if (line.size() >= 51) {
					valence[i] = lexical_cast<Index>(strip(line.substr(48, 3)));
				}
				if (line.size() >= 63) {
					aam[i] = lexical_cast<Index>(strip(line.substr(60, 3)));
				}
				if (line.size() >= 66) {
					inversion[i] = lexical_cast<Index>(strip(line.substr(63, 3)));
				}
				if (symbol == "D") {
					ct0.back() = 1;
					iso[i] = 2;
				}
				else if (symbol == "T") {
					ct0.back() = 3;
					iso[i] = 2;
				}
			}
		}

		for (size_t i = 0; i < nBonds; ++i) {
			if (getline(stream, line)) {
				Index a1 = lexical_cast<Index>(strip(line.substr(0, 3)))-1;
				Index a2 = lexical_cast<Index>(strip(line.substr(3, 3)))-1;
				if (a1 >= nAtoms || a2 >= nAtoms) {
					throw error("The atom index is out of range");
				}
				Index o = lexical_cast<Index>(strip(line.substr(6, 3)));
				switch (o) {
				case 1:
					o = Index(Bond::SINGLE_BOND);
					break;

				case 2:
					o = Index(Bond::DOUBLE_BOND);
					break;

				case 3:
					o = Index(Bond::TRIPLE_BOND);
					break;

				case 4:
					o = Index(Bond::PARTIAL_BOND);
					break;
				}
				ct0.push_back(o);
				ct0.push_back(a1);
				ct0.push_back(a2);
				order[i] = o;

				if (line.size() >= 12) {
					stereoB[i] = lexical_cast<Index>(strip(line.substr(9, 3)));
				}
				if (line.size() >= 21) {
					center[i] = lexical_cast<Index>(strip(line.substr(18, 3)));
				}
			}
		}

		if (stream) {
			while (getline(stream, line)) { 
				boost::algorithm::trim(line);
				if (line == "M  END") {
					break;
				}
				else if (boost::algorithm::starts_with(line, "M  CHG")) {
					size_t n = lexical_cast<size_t>(boost::algorithm::trim_copy(line.substr(6, 3)));
					for (size_t i = 0, offset=10; i < n; ++i, offset += 8) {
						Index ia = lexical_cast<Index>(boost::algorithm::trim_copy(line.substr(offset, 4)))-1;
						Integer chg = lexical_cast<Integer>(boost::algorithm::trim_copy(line.substr(offset+4, 4)));
						if (ia < charge.size()) {
							charge[ia] = chg;
						}
						else {
							throw BadLexicalCast(line, typeid(Atom));
						}
					}
				}
				else if (boost::algorithm::starts_with(line, "M  RAD")) {
					size_t n = lexical_cast<size_t>(boost::algorithm::trim_copy(line.substr(6, 3)));
					for (size_t i = 0, offset=10; i < n; ++i, offset += 8) {
						Index ia = lexical_cast<Index>(boost::algorithm::trim_copy(line.substr(offset, 4)))-1;
						Integer radI = lexical_cast<Integer>(boost::algorithm::trim_copy(line.substr(offset+4, 4)));
						if (ia < rad.size()) {
							rad[ia] = radI;
						}
						else {
							throw BadLexicalCast(line, typeid(Atom));
						}
					}				
				}
				else if (boost::algorithm::starts_with(line, "M  ISO")) {
					size_t n = lexical_cast<size_t>(boost::algorithm::trim_copy(line.substr(6, 3)));
					for (size_t i = 0, offset=10; i < n; ++i, offset += 8) {
						Index ia = lexical_cast<Index>(boost::algorithm::trim_copy(line.substr(offset, 4)))-1;
						Integer isoI = lexical_cast<Integer>(boost::algorithm::trim_copy(line.substr(offset+4, 4)));
						if (ia < iso.size()) {
							iso[ia] = isoI;
						}
						else {
							throw BadLexicalCast(line, typeid(Atom));
						}
					}				
				}

			}

			assert (ct0.size() == 2 + nAtoms + 3*nBonds);

			result[MOL_INFO_CT] = SimObjPtr(new PropertyConnectTable(&ct0[0]));
			result[MOL_INFO_ATOM_COORD_2D] = SimObjPtr(new PropertyXYZ(coords, false));
			result[MOL_INFO_ATOM_COORD_3D] = SimObjPtr(new PropertyXYZ(coords, false));

			if (!isZero(iso)) {
				result[MOL_INFO_ATOM_ISO] = SimObjPtr(new PropertyIndex(iso, true));
			}
			if (!isZero(charge)) {
				result[MOL_INFO_ATOM_CHARGE] = SimObjPtr(new PropertyInteger(charge, true));
			}
			if (!isZero(rad)) {
				result[MOL_INFO_ATOM_RADICAL] = SimObjPtr(new PropertyIndex(rad, true));
			}
			if (!isZero(stereoA)) {
				// result[MOL_INFO_ATOM_STEREO] = SimObjPtr(new PropertyIndex(stereoA, true));
			}
			if (!isZero(implicitH)) {
				for (size_t i = 0; i < implicitH.size(); ++i) {
					implicitH[i] -= 1;
				}
				result[MOL_INFO_ATOM_IMPLICIT_H] = SimObjPtr(new PropertyIndex(implicitH, true));
			}
			if (!isZero(valence)) {
				// result[MOL_INFO_ATOM_VALENCE] = SimObjPtr(new PropertyIndex(valence, true));
			}
			if (!isZero(aam)) {
				result[RXN_INFO_ATOM_AAM] = SimObjPtr(new PropertyIndex(aam, 0));
			}
			if (!isZero(inversion)) {
				result[RXN_INFO_ATOM_INVERSION] = SimObjPtr(new PropertyIndex(inversion, true));
			}

			if (!isZero(stereoB)) {
				result[MOL_INFO_BOND_STEREO] = SimObjPtr(new PropertyIndex(stereoB, true));
			}

			if (!isZero(center)) {
				result[RXN_INFO_BOND_CENTER] = SimObjPtr(new PropertyIndex(center, true));
			}

			return true;
		}
		else {
			return false;
		}
	};


	bool MDLMolWriter::write(std::ostream& stream, const PhaseSpace& data) const {

		stream << "\n    Kuai Mol\n\n";

		if (data.find(MOL_INFO_CT) != data.end()) {
			const ConnectTable& ct = static_cast<const PropertyConnectTable&>(*(data.find(MOL_INFO_CT)->second));
			std::vector<Index> zero(std::max(ct.countAtoms(), ct.countBonds()), 0);
			std::vector<XYZ> zeroXYZ(ct.countAtoms(), XYZ(0, 0, 0));

			const XYZ* pCoord = &zeroXYZ[0];
			if (data.find(MOL_INFO_ATOM_COORD_3D) != data.end()) {
				pCoord = static_cast<const PropertyXYZ&>(*(data.find(MOL_INFO_ATOM_COORD_3D)->second)).c_ptr();
			}
			else if (data.find(MOL_INFO_ATOM_COORD_2D) != data.end()) {
				pCoord = static_cast<const PropertyXYZ&>(*(data.find(MOL_INFO_ATOM_COORD_2D)->second)).c_ptr();
			}

			const Index* pOrder = NULL;
			if (data.find(MOL_INFO_BOND_ORDER) != data.end()) {
				pOrder = static_cast<const PropertyIndex&>(*(data.find(MOL_INFO_BOND_ORDER)->second)).c_ptr();
			}

			char buf[128];
			
			sprintf(buf, "%3d%3d  0  0  0  0  0  0  0  0999 V2000\n", ct.countAtoms(), ct.countBonds());
			stream << buf;
			for (size_t i = 0; i < ct.countAtoms(); ++i) {
				sprintf(buf, "%10.4f%10.4f%10.4f %-2s  0  0  0  0  0  0  0  0  0  0  0  0\n",
					pCoord[i].x, pCoord[i].y, pCoord[i].z, Element::getElement(ct.getAtom(i).type).symbol.c_str());
				stream << buf;
			}
			for (size_t i = 0; i < ct.countBonds(); ++i) {
				const Bond& bondI = ct.getBond(i);
				Index a1 = ct.indexOf(ct.getAtom(bondI, 0));
				Index a2 = ct.indexOf(ct.getAtom(bondI, 1));
				
				Index o;
				if (pOrder) {
					o = pOrder[i];
				}
				else {
					o = bondI.order;
				}
				switch (o) {
				case Bond::SINGLE_BOND:
					o=1;
					break;
				case Bond::DOUBLE_BOND:
					o=2;
					break;
				case Bond::TRIPLE_BOND:
					o=3;
					break;
				case Bond::PARTIAL_BOND:
					o=4;
					break;
				default:
					o=0;
					break;
				}

				sprintf(buf, "%3d%3d%3d  0  0  0  0\n",
					a1+1, a2+1, o);
				stream << buf;
			}
			stream << "M  END\n";
		}
	
		return true;
	}

	
	namespace {
	
		static const Char SZ_DEFAULT_KUAIMOL_ENCRYPT_CODE[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz()=";
	
		inline Index width(Index nAtoms, Index nBonds) {
			Index v = std::max(nAtoms, nBonds);
			Index width = 0;
			while (v > 0) {
				v /= 64;
				++width; 
			}
			return width;
		}
		inline void dumpIndex(std::ostream& stream, Index id, Index width, const Char code[]) {
			for (Index i = 0; i < width; ++i) {
				stream << code[id % 64];
				id /= 64;
			}
			assert (id == 0);
		}
		inline void dumpRealNumber(std::ostream& stream, RealNumber v, const Char code[]) {
			assert (-8388 < v && v < 8388);
			boost::uint32_t vv = boost::uint32_t((v + 8388) * 1000);
			for (Index i = 0; i < 4; ++i) {
				stream << code[vv % 64];
				vv /= 64;
			}
			assert (vv == 0);
		}
		inline RealNumber loadRealNumber(std::istream& stream, const Char code[]) {
			boost::uint32_t result = 0;
			boost::uint32_t shift = 1;
			for (Index i = 0; i < 4; ++i) {
				std::istream::int_type c = stream.get();
				if (c < 0  || 128 <= c) {
					stream.setstate(std::ios_base::failbit);
					return 0;
				} 
				else {
					result += shift * code[c];
					shift *= 64;
				}
			}
			return RealNumber(result) / 1000 - 8388;
		}
		inline Index loadIndex(std::istream& stream, Index width, const Char code[]) {
			Index result = 0;
			Index shift = 1;
			for (Index i = 0; i < width; ++i) {
				std::istream::int_type c = stream.get();
				if (c < 0  || 128 <= c) {
					stream.setstate(std::ios_base::failbit);
					return Index(-1);
				} 
				else {
					result += code[c] * shift;
					shift *= 64;
				}
			}
			return result;
		}

		class KuaiMolIOManager {
		public:
			typedef bool (*PropertySaver)(std::ostream&, const SimObject& obj, Index width, const Char szCoder[]);
			typedef bool (*PropertyLoader)(std::istream&, SimObject& obj, Index width, const Char szCoder[]);

		private:
			KuaiMolIOManager() {
				addHandle(MOL_INFO_ATOM_COORD_2D, saveCoord2D, loadCoord2D);
				addHandle(MOL_INFO_ATOM_COORD_3D, saveCoord3D, loadCoord3D);
				addHandle(MOL_INFO_BOND_ORDER, saveIndexArray1, loadIndexArray1);
				addHandle(MOL_INFO_BOND_RANK, saveIndexArray, loadIndexArray);

				// addHandle(MOL_INFO_ATOM_CHARGE);
				// addHandle(MOL_INFO_ATOM_ISO);
				// addHandle(MOL_INFO_ATOM_RADICAL);
				// addHandle(MOL_INFO_ATOM_STEREO);
				// addHandle(MOL_INFO_ATOM_IMPLICIT_H);
				// addHandle(MOL_INFO_ATOM_VALENCE);
				
				// addHandle(MOL_INFO_BOND_STEREO);
				
				// addHandle(RXN_INFO_ATOM_AAM);
				// addHandle(RXN_INFO_ATOM_INVERSION);
				// addHandle(RXN_INFO_BOND_CENTER);
			}

		public:
			static KuaiMolIOManager& getInstance() {
				static KuaiMolIOManager singleton;
				return singleton;
			}
			bool hasHandle(const String& prop) const {
				return prop2type.find(prop) != prop2type.end();
			};
			bool hasHandle(Index id) const {
				return id < type2prop.size();
			};

			bool addHandle(const String& prop, PropertySaver saver, PropertyLoader loader) {
				if (hasHandle(prop)) {
					return false;
				}
				else {
					prop2type[prop] = type2prop.size();
					type2prop.push_back(prop);
					savers.push_back(saver);
					loaders.push_back(loader);
					return true;
				}
			}

			bool save(std::ostream& stream, PhaseSpace::const_iterator it, Index width, const Char szCoder[]) const {
				assert (hasHandle(it->first));
				Index i = prop2type.find(it->first)->second;
				dumpIndex(stream, i, 2, szCoder);
				return savers[i](stream, *(it->second), width, szCoder);
			}

			bool load(std::istream& stream, Index id, PhaseSpace& space, Index width, const Char szCoder[]) const {
				assert (hasHandle(id));
				SimObjPtr ptr = space[type2prop[id]];
				
				if (loaders[id](stream, *ptr, width, szCoder)) {
					return true;
				}
				else {
					return false;
				}
			}

		private:
			std::vector<String>   type2prop;
			std::map<String, Index> prop2type;
			std::vector<PropertySaver> savers;
			std::vector<PropertyLoader> loaders;

		private:
			static bool saveCoord2D(std::ostream& stream, const SimObject& obj, Index width, const Char szCoder[]) {
				const PropertyXYZ& xyz = static_cast<const PropertyXYZ&>(obj);
				size_t n = xyz.size();
				const XYZ* p = xyz.c_ptr();
				for (size_t i = 0; i < n && stream.good(); ++i) {
					dumpRealNumber(stream, p[i].x, szCoder);
					dumpRealNumber(stream, p[i].y, szCoder);
				}
				return stream.good();
			}
			static bool loadCoord2D(std::istream& stream, SimObject& obj, Index width, const Char szCoder[]) {
				PropertyXYZ& xyz = static_cast<PropertyXYZ&>(obj);
				size_t n = xyz.size();
				XYZ* p = xyz.x_ptr();
				for (size_t i = 0; i < n && stream.good(); ++i) {
					p[i].x = loadRealNumber(stream, szCoder);
					p[i].y = loadRealNumber(stream, szCoder);
					p[i].z = 0;
				}
				return stream.good();
			}
			static bool saveCoord3D(std::ostream& stream, const SimObject& obj, Index width, const Char szCoder[]) {
				const PropertyXYZ& xyz = static_cast<const PropertyXYZ&>(obj);
				size_t n = xyz.size();
				const XYZ* p = xyz.c_ptr();
				for (size_t i = 0; i < n && stream.good(); ++i) {
					dumpRealNumber(stream, p[i].x, szCoder);
					dumpRealNumber(stream, p[i].y, szCoder);
					dumpRealNumber(stream, p[i].z, szCoder);
				}
				return stream.good();
			}
			static bool loadCoord3D(std::istream& stream, SimObject& obj, Index width, const Char szCoder[]) {
				PropertyXYZ& xyz = static_cast<PropertyXYZ&>(obj);
				size_t n = xyz.size();
				XYZ* p = xyz.x_ptr();
				for (size_t i = 0; i < n && stream.good(); ++i) {
					p[i].x = loadRealNumber(stream, szCoder);
					p[i].y = loadRealNumber(stream, szCoder);
					p[i].z = loadRealNumber(stream, szCoder);
				}
				return stream.good();
			}

			static bool saveIndexArray(std::ostream& stream, const SimObject& obj, Index width, const Char szCoder[]) {
				const PropertyIndex& v = static_cast<const PropertyIndex&>(obj);
				size_t n = v.size();
				const Index* p = v.c_ptr();
				for (size_t i = 0; i < n && stream.good(); ++i) {
					dumpIndex(stream, p[i], width, szCoder);
				}
				return stream.good();
			}
			static bool loadIndexArray(std::istream& stream, SimObject& obj, Index width, const Char szCoder[]) {
				PropertyIndex& v = static_cast<PropertyIndex&>(obj);
				size_t n = v.size();
				Index* p = v.x_ptr();
				for (size_t i = 0; i < n && stream.good(); ++i) {
					p[i] = loadIndex(stream, width, szCoder);
				}
				return stream.good();
			}

			static bool saveIndexArray1(std::ostream& stream, const SimObject& obj, Index width, const Char szCoder[]) {
				const PropertyIndex& v = static_cast<const PropertyIndex&>(obj);
				size_t n = v.size();
				const Index* p = v.c_ptr();
				for (size_t i = 0; i < n && stream.good(); ++i) {
					dumpIndex(stream, p[i], 1, szCoder);
				}
				return stream.good();
			}
			static bool loadIndexArray1(std::istream& stream, SimObject& obj, Index width, const Char szCoder[]) {
				PropertyIndex& v = static_cast<PropertyIndex&>(obj);
				size_t n = v.size();
				Index* p = v.x_ptr();
				for (size_t i = 0; i < n && stream.good(); ++i) {
					p[i] = loadIndex(stream, 1, szCoder);
				}
				return stream.good();
			}
		};
	}

	KuaiMolReader::KuaiMolReader() {
		code.resize(128, 0);
		for (Index i = 0; i < 64; ++i) {
			code[SZ_DEFAULT_KUAIMOL_ENCRYPT_CODE[i]] = Char(i);
		}
		TAG_OF_END_OF_BLOCK = SZ_DEFAULT_KUAIMOL_ENCRYPT_CODE[64];
	}
	KuaiMolReader::KuaiMolReader(const Char szCoder[]) {
		code.resize(128, 0);
		for (Index i = 0; i < 64; ++i) {
			code[szCoder[i]] = Char(i);
		}
		TAG_OF_END_OF_BLOCK = szCoder[64];
	}
	bool KuaiMolReader::update(std::istream& stream, PhaseSpace& space) const {
		Char buffer[64];
		if (!stream.read(buffer, 4) || memcmp(buffer, "[M00", 4) != 0) {
			return false;
		}
		PropertyConnectTable* pCT= readConnectTable(stream);
		if (pCT == NULL || !stream) {
			return false;
		}
		space[MOL_INFO_CT] = SimObjPtr(pCT);
		Index w = width(pCT->c_ct().countAtoms(), pCT->c_ct().countBonds());

		stream.get();	// Ignore the TAG_OF_END_OF_BLOCK;
		for (;;) {
			if (stream.read(buffer, 2)) {
				if (memcmp(buffer, "M]", 2) == 0) {
					break;
				}
				else {
					Index blockID = code[buffer[0]] + code[buffer[1]]*64;
					if (KuaiMolIOManager::getInstance().hasHandle(blockID)) {
						KuaiMolIOManager::getInstance().load(stream, blockID, space, w, &code[0]);
						stream.get();
					}
					else {
						String temp;
						std::getline(stream, temp, TAG_OF_END_OF_BLOCK);
					}
				}
			}
			else {
				break;
			}
		}
		
		return stream.good();
	}

	PropertyConnectTable* KuaiMolReader::readConnectTable(std::istream& stream) const {
		Index nAtoms = loadIndex(stream, 3, &code[0]);
		Index nBonds = loadIndex(stream, 3, &code[0]);
		if (stream.good()) {
			std::vector<Index> v; v.reserve(2 + nAtoms + 3 * nBonds);
			v.push_back(nAtoms), v.push_back(nBonds);
			Index w = width(nAtoms, nBonds);
			Index nHAtTail = loadIndex(stream, w, &code[0]);
			if (nHAtTail <= nAtoms) {
				{
					Index i;
					for (i= 0; i < nHAtTail; ++i) {
						v.push_back(6);
					}
					for (; i < nAtoms; ++i) {
						v.push_back(1);
					}
				}

				for (;;) {
					Index i = loadIndex(stream, w, &code[0]);
					if (i >= nAtoms) {
						break;
					}
					else {
						assert (i < nHAtTail);
						v[i+2] = loadIndex(stream, 2, &code[0]);
					}
				}

				for (Index i = 0; i < nBonds; ++i) {
					Index a1 = loadIndex(stream, w, &code[0]);
					Index a2 = loadIndex(stream, w, &code[0]);
					Index order = loadIndex(stream, 1, &code[0]) / 8;

					if (a1 != a2 && a1 < nAtoms && a2 < nAtoms && order <= Bond::TRIPLE_BOND) {
						v.push_back(order);
						v.push_back(a1);
						v.push_back(a2);
					}
					else {
						stream.setstate(std::ios_base::failbit);
						break;
					}
				}

				if (stream.good()) {
					assert (v.size() == 2 + nAtoms + 3 * nBonds);
					return new PropertyConnectTable(&v[0]);
				}
			}
		}
		return NULL;
	};

	KuaiMolWriter::KuaiMolWriter() {
		code = SZ_DEFAULT_KUAIMOL_ENCRYPT_CODE;
	}

	KuaiMolWriter::KuaiMolWriter(const Char szCode[]) 
		: code(szCode)
	{
		assert (code.size() == 65);
	}


	bool KuaiMolWriter::write(std::ostream& stream, const PhaseSpace& data) const {
		const ConnectTable& ct = data.ct();
		stream << "[M00";
		Index width = writeConnectTable(stream, ct);
		stream << code[64];
		for (PhaseSpace::const_iterator 
			it = data.begin(); it != data.end(); ++it)
		{
			if (KuaiMolIOManager::getInstance().hasHandle(it->first)) {
				if (!KuaiMolIOManager::getInstance().save(stream, it, width, code.c_str())) {
					break;
				}
				stream << code[64];
			}
		}
		return (stream << "M]").good();
	}
	
	Index KuaiMolWriter::writeConnectTable(std::ostream& stream, const ConnectTable& ct) const {
		Index nAtoms = ct.countAtoms();
		Index ait = 0;
		while (nAtoms > 0) {
			stream << code[nAtoms % 64];
			++ait;
			nAtoms /= 64;
		}
		for (Index i = ait; i < 3; ++i) {
			stream << code[0];
		}
		nAtoms = ct.countAtoms();

		Index nBonds = ct.countBonds();
		Index bit = 0;
		while (nBonds > 0) {
			stream << code[nBonds % 64];
			++bit;
			nBonds /= 64;
		}
		for (Index i = bit; i < 3; ++i) {
			stream << code[0];
		}
		nBonds = ct.countBonds();

		Index width = std::max(ait, bit);

		Index nHAtTail = 0;
		for (nHAtTail = nAtoms+1; nHAtTail > 0;) {
			if (ct.getAtom(--nHAtTail).type != 1) {
				break;
			}			
		}
		dumpIndex(stream, nHAtTail, width, code.c_str());

		for (Index i = 0; i < nHAtTail; ++i) {
			const Atom& atomI = ct.getAtom(i);
			if (atomI.type != 6) {
				assert (atomI.type < 64*64);
				dumpIndex(stream, i, width, code.c_str());
				dumpIndex(stream, atomI.type, 2, code.c_str());
			}
		}
		dumpIndex(stream, nAtoms, width, code.c_str());

		for (Index i = 0; i < nBonds; ++i) {
			const Bond& bondI = ct.getBond(i);
			dumpIndex(stream, ct.indexOf(ct.getAtom(bondI, 0)), width, code.c_str());
			dumpIndex(stream, ct.indexOf(ct.getAtom(bondI, 1)), width, code.c_str());
			dumpIndex(stream, bondI.order * 8 + rand() % 8, 1, code.c_str());
		}

		return width;
	}

	MolReaderManager::MolReaderManager() {
		(*this)[".mol"] = (*this)[".mdl"] = boost::shared_ptr<SpaceReader>(new MDLMolReader);
		(*this)[".kmf"] = boost::shared_ptr<SpaceReader>(new KuaiMolReader);
		// (*this)[".gro"] = boost::shared_ptr<SpaceReader>(new GromacsGROReader);
		// (*this)[".xyz"] = boost::shared_ptr<SpaceReader>(new XYZReader);
	}

	MolReaderManager& MolReaderManager::getInstance() {
		static MolReaderManager singleton;
		return singleton;
	};

	MolWriterManager::MolWriterManager() {
		(*this)[".mol"] = (*this)[".mdl"] = boost::shared_ptr<SpaceWriter>(new MDLMolWriter);
		(*this)[".kmf"] = boost::shared_ptr<SpaceWriter>(new KuaiMolWriter);
	}

	MolWriterManager& MolWriterManager::getInstance() {
		static MolWriterManager singleton;
		return singleton;
	};

	const char* getDefaultKuaiMolKey() {
		return SZ_DEFAULT_KUAIMOL_ENCRYPT_CODE;	
	};
}
